A little background on work done on the tel may prove helpful.
Information of the mapping of the tel can be found at http://facweb1.redlands.edu/fac/jim_bentley/talks/elFaraSMap.ppt
Information of the surface survey can be found at http://facweb1.redlands.edu/fac/jim_bentley/talks/elFaraSSurvey.ppt
Get the sherds data from the web.
circles = read.csv("http://facweb1.redlands.edu/facultyfolder/jim_bentley/downloads/math312/circles9812.csv")
names(circles)
## [1] "EFT" "EM" "NFT" "NM" "HFT" "HM"
## [7] "SESSION" "ESTEM" "ESTNM" "ESTHM" "THEOCODE" "DATE"
## [13] "AREA" "COLUNIT" "RADIUS" "MB" "LB" "IAI"
## [19] "IAII" "IAIII" "PERS" "HELL" "ROM" "BYZ"
## [25] "MA" "OTT" "UNIDENT" "TOTAL" "PH" "HR"
## [31] "PR" "NOTE"
dim(circles)
## [1] 70 32
From the above we can see that there are 32 variables and 70 observations. Two of the variables contain information on the number of iron age II sherds found in each circle and the easting of the center of the circle in meters.
circles$IAII
## [1] 0 3 0 0 0 2 0 1 6 0 0 2 0 0 0 0 0 1 0 1 3 2 1
## [24] 2 2 0 0 0 0 0 5 5 2 0 0 0 0 0 0 0 16 8 0 7 0 2
## [47] 0 8 2 0 0 1 0 7 5 3 5 1 0 1 0 7 4 12 3 0 0 1 1
## [70] 3
circles$ESTEM
## [1] NA NA 100703.6 100704.3 100693.1 100682.7 100698.8
## [8] 100678.4 100669.1 100686.2 100668.4 100692.6 100689.6 100678.1
## [15] 100683.0 100671.8 100670.4 100695.5 100688.4 100660.5 100648.8
## [22] 100672.0 100675.7 100663.4 100639.2 100641.6 100625.3 100624.9
## [29] 100630.4 100635.5 100628.1 100613.9 100609.7 100594.4 100598.3
## [36] 100654.4 100643.1 100647.0 100630.3 100634.8 100624.3 100675.1
## [43] 100667.0 100681.8 100676.2 100665.8 100652.0 100655.1 100688.0
## [50] 100656.2 100696.5 100702.3 100674.6 100636.6 100640.9 100597.7
## [57] 100606.6 100594.6 100590.3 100598.5 100585.5 100596.0 100610.9
## [64] 100627.0 100637.2 100642.0 100608.6 100628.5 100641.7 100618.0
length(circles$IAII)
## [1] 70
length(circles$ESTEM)
## [1] 70
The first two rows are trash. We can remove them.
circles=circles[-c(1:2),]
dim(circles)
## [1] 68 32
We can let our \(X_i\) represent the Iron Age II sherds.
x = circles$IAII
table(x)
## x
## 0 1 2 3 4 5 6 7 8 12 16
## 34 9 8 4 1 4 1 3 2 1 1
The circles with no sherds are meaningless as they do not affect the sum. We remove them.
x.no0 = x[x != 0]
length(x.no0)
## [1] 34
x.no0
## [1] 2 1 6 2 1 1 3 2 1 2 2 5 5 2 16 8 7 2 8 2 1 7 5
## [24] 3 5 1 1 7 4 12 3 1 1 3
We can now compute \(\Lambda\) and the associated p-values.
(n <- length(x))
## [1] 68
(df <- n-1)
## [1] 67
(Lambda <- 2*t(x.no0)%*%log(x.no0/mean(x)))
## [,1]
## [1,] 267.0377
1-pchisq(Lambda,df)
## [,1]
## [1,] 0
(test.stat <- (n-1)*sqrt(var(x))/mean(x))
## [1] 107.3201
1-pchisq(test.stat,df)
## [1] 0.001288988
A variogram helps us see how the sherd counts vary across the tel.
dists=dist(circles[,c("ESTEM","ESTNM")])
summary(dists)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.113 45.916 74.171 77.921 105.964 188.652
breaks=seq(0,190,l=11)
vl = variog(coords=circles[,c("ESTEM","ESTNM")],data=circles$IAII,breaks=breaks)
## variog: computing omnidirectional variogram
vl.summary = cbind(c(1:10),vl$v,vl$n)
colnames(vl.summary) = c("lag","semi-variance","# of pairs")
vl.summary
## lag semi-variance # of pairs
## [1,] 1 5.692661 109
## [2,] 2 10.593220 295
## [3,] 3 9.000000 373
## [4,] 4 9.597744 399
## [5,] 5 10.277778 351
## [6,] 6 8.878333 300
## [7,] 7 9.400000 210
## [8,] 8 10.964029 139
## [9,] 9 14.388889 81
## [10,] 10 7.500000 21
plot(vl,type="b")